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ABSTRACT 

We consider travel time tomography problems involving detection of high contrast, discrete 
high velocity structures. This results in a discrete nonlinear inverse problem, for which tra- 
ditional grid-based models and iterative linearized least-squares reconstruction algorithms 
are not suitable. This is because travel paths change significantly near the high contrast 
velocity structure, making it more difficult to inversely calculate the travel path and in- 
fer the velocity along the path. We propose a model-based approach to describe the high 
velocity structure using pre-defined elementary objects. Compared to a grid-based model, 
our approach has complexity that increases as a function of the number of objects, rather 
than increasing with the number of cells (usually very large). A new reconstruction algo- 
rithm is developed that provides estimates of the probability that a high velocity structure 
appears at any point in the region of interest. Simulation results show that our method can 
efficiently sample the model parameter space, and we map the model parameters into the 
high velocity structures in spatial domain to generate a "probability map" , which represent 
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the appearance of the high velocity structure in different regions. We show the probabihty 
map not only gives the highest probability to the optimal solution, but also includes other 
possible models as well. 
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INTRODUCTION 

Travel time tomography aims to reconstruct an interior velocity model by using the mea- 
sured first-arrival times between transmitters and receivers. The velocity model captures 
the physical properties of the region where the signal transmission occurs. The transmit- 
ted signal can be a seismic wave, an acoustic sound wave, and even a fluid pressure wave. 
This technique is important to characterize large scale elastic media, in applications such 



as seismic geophysical exploration (Pratt, 1999) and acoustic tomography (for example 



in the atmosphere's surface (Wilson et al., 2001) or in the ocean layers (Munk and Wun 



sch, 1979[ )). However, different from many other transmission tomographic reconstruction 



problems (X-Ray CT, Positron emission tomography CT) where the straight-line trajec- 
tory assumption is commonly used in the problems we consider here the travel path may 
bend severely if the travel velocity variation is very high. This phenomenon happens in all 
kinds of waveform inversion applications, with one of the most studied cases being acoustic 
tomography. Instead of using the straight line trajectory assumption, travel paths can be 
characterized by Fermat's principle: the travel path observed is the path traversed in the 



shortest time (Berryman, 1989), where the time can be predicted as a path integral of the 



slowness function (reciprocal of velocity) integrated along the travel path. 

In many geophysical applications it is very common to have heterogeneous structures 
where the velocity contrast is high, e.g., a factor of two difference in velocity between 
different areas. Example of scenarios where this situation can be encountered arise in many 
applications: using seismic waves to find permeable fracture zones in ground-water flow 



characterization (National Research Council (US). Committee on Fracture Characterization 



and Fluid Flow, 1996), monitoring the water/oil saturation in vegetable oil bio-remediation 
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projects (Lane Jr et al., 2004), etc. Different from the case where the velocity distribution 



only has small variations (|Wilson et al., 2001), the travel path in high velocity contrast 



media not only bends severely but is almost dominated by these high velocity structures 
that form high transmissive channels. Because the straight line travel path approximation 
is not valid, the reconstruction of the velocity model becomes a nonlinear inverse problem. 
Note that conventional reconstruction methods are based on iterative linearized algorithms 
to approximate the travel path. Thus, given that in many practical situations we only 
have sparse measured travel time data, these methods suffer from problems due to low ray- 
coverage and severe path bending in the high contrast velocity medium. For a comprehensive 



review on this subject, we refer the reader to the overview by Berryman (Berryman, 1991a). 



Our initial motivation for this work comes from the flow permeability characterization 



in a fractured reservoir (Vasco and Datta-Gupta, 1997). One of the most widely used 



enhanced oil recovery (EOR) techniques, water-flooding, involves injecting water in a con- 
trolled manner in order to provide pressure support that can slowly sweep the oil into the 



production wells (Welge, 1952). During this process, the permeability (measurement of 
ability to transmit fluid) of open fractures can be orders of magnitude higher than that of 
surrounding tight rocks, providing fast pathways for the flow. Thus, travel time through 
a fracture (which could be modeled as a line in 2D or a plane in 3i?) is much faster than 
through surrounding rocks. The flow properties of the reservoir are dominated by these 
highly 'transmissive' fracture structures. If a fracture is close to both an injection well and 
a production well, most water will flow directly through the fracture and will fail to displace 
oil in nearby areas. This is a phenomenon known as "water cycling", which significantly 
reduces oil recovery efficiency. Thus, understanding the locations of fractures is critical in 
flow characterization and enhancement of oil recovery efficiency. The major challenge to 



achieve this understanding is that travel time measurements are hmited by the bore-hole lo- 
cations, which makes it difficult to reconstruct high resolution images of fracture locations. 
The high permeability contrast property between open fractures and surrounding rocks also 
increases the difficulty of reconstructing accurate velocity maps. 

A similar situation also arises in bio-remediation problems, where the goal is to map 
the hydraulic conductivity and predict the ground water flow in the aquifer. Bore-hole core 
samples and pollutant concentration data are available only from the pumping wells, which 
are very sparsely located in the fleld. Another interesting application is network delay 



tomography (Tsang et al., 2003), where the main goal is to identify the nodes where con- 
gestion occurs. Instead of using internal measurements at many nodes, which require high 
bandwidth resources, network tomography measures the end-to-end delay in the network 
and estimates the delay in each link. When congestion happens, the delay time may be 
very different from the typical values, which would make this a 'high-contrast' tomography 
problem. 

Many different models have been developed for problems that involve a structured ve- 



locity distribution. Grid-based methods (Berryman, 1991b), which divide the volume into 



small cells and assign a constant velocity to each cell, are widely used because of the freedom 
to represent the structure in any degree of detail by increasing the number of cells. As an 



alternative, object-based methods (Lane Jr et al., 2004) use predeflned objects to represent 



the velocity structure. Compared to grid-based methods, if the pre-deflned shape of those 
objects is chosen wisely, object-based methods have the advantage of parameterizing the 
spatial distribution of velocity with a small number of objects instead of a large number of 
cells. 
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This paper focuses on finding high velocity structures (HVS) in a relatively slow homo- 
geneous background. Specifically we consider the case where the travel velocity can only 
take one of several possible discrete values. We call this a 'discrete' travel time tomography 
problem. For example, a reservoir can usually be modeled as a layered structure, where 
each layer has quite different hydraulic conductivity. 

Only a limited amount of research has addressed the high contrast velocity problem in 



travel time tomography. Bai and Greenhalgh (Bai and Greenhalgh, 2005) proposed to add 



irregular nodes on the boundary of cells to improve the stability when determining nonlinear 



ray travel paths. Berryman (Berryman, 1990) used Fermat's principle on the velocity model 



to handle the case where high contrast velocity alters the travel path severely. Both of these 
approaches are grid-based methods and use iterative linearized algorithms to perform the 
tomographic reconstruction. Because the velocity contrast is so high, grid-based models 
require a very fine grid to represent the velocity distribution at the boundaries between 
areas of different velocity, and the corresponding iterative linearized inversion algorithm 
will often fail to converge and determine the actual travel path. A finer grid implies we 
need to estimate more unknown values (the number of parameters to estimate grows linearly 
with the number of cells), some of which may not even be covered by any travel path. For 
these areas, we cannot determine the travel velocity because no travel path passes through 
the corresponding cells. This is the well known "lack of ray coverage" phenomenon. 

As an alternative, in this paper, we use object-based models to represent the HVS. The 
objects are chosen from pre-defined fundamental convex shapes, so that the geometrical 
shape of HVS can be represented by a combination of these objects. This approach has 
two main advantages. First, the problem of approximating an arbitrary shape by multiple 
convex objects is well understood. Moreover, we can incorporate the prior information 
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about the HVS into object-based models. For instance, in the fracture characterization 
problem we know that the geometrical shape of fractures can be well approximated by lines 
in 2D models or planes in 31? models. Thus, in these cases, we can define the shape of our 
fundamental convex object as a "line" or a "plane", respectively. With this formulation, 
only small number of objects is needed to well approximate the shape of HVS, which reduces 
the dimensionality and uncertainty of the reconstruction problem. Second, the travel path 
tracking procedure can be simplified by only considering the shortest path between objects 
instead of all cells in spatial domain. Because the elementary objects are convex, the 
travel path is piecewise linear and the computational complexity for finding the shortest 
path between objects is reduced. Compared to the methods we mentioned before, our 
approach reduces the number of unknown variables (which is equal to the number of object 
parameters) and avoids the "lack of path coverage" problem arising in grid-based models, 
which leads to a more stable solution to the problem. 

To estimate the velocity model from the measured travel time, we need to solve a 
nonlinear inverse problem. The major issue in all inverse problems is that the solution, in 
this case the estimated velocity model, may not be unique for the given measurements. One 
popular approach to handle the non-uniqueness issue is to apply regularizations to favor 



certain properties in the model (Engl et al., 1996). The regularization methods can be 



viewed as model selection: they will lead to solutions that balance data-fitting and model- 
penalty. However, it is not easy to select the weight for model-penalty and this usually 



requires cross validation in order to avoid over-fitting (Ng, 1997). 



An alternative approach is to estimate the probability distribution for the model space 



according to the data- fitting (Tarantola, 1987). This gives a full description of the relative 



probabilities of the different models, so that we can explicitly consider all likely solutions. 



However, generating samples and estimating the probability density is very challenging and 



time consuming for a high dimensional model space (MacKay, 2003). 



In this paper, we choose the Bayesian approach and focus on estimating the probability 
density in the model parameter space. Due to the large dimension of the model parameter 
space, we develop an accelerated random walk algorithm to explore the model space. Based 
on the Hamiltonian Monte Carlo method (HMC), we add an additional friction term to 
the generation of sampling points along the simulated particle moving trajectory through 
the model parameter space. Our algorithm will sample more frequently locally in low 
error regions, and re-sampling will be used based on the HVS property to approximate 
the structure of the probability distribution. The results are presented as a probability 
map showing where the high velocity objects are more likely to appear in the underlying 
structure of the system. 

To the best of our knowledge, we are the first to use an object-based approach to solve 
a high contrast, discrete velocity tomography problem. Our proposed algorithm uses the 
HVS properties to simplify the path finding step, which makes the HMC sampling process 
more efficient. Furthermore, we exploit the monotonicity of the travel time as a function of 
high velocity object size to approximate the probability distribution in the space of model 
parameters. The rest of the paper is organized as follows: We define the object based 
model to represent the structure in the velocity model in the next section, and introduce the 
forward step and define the mathematical formulation for the travel path finding problem. 
Then we give an overview of our proposed algorithm for probability density estimation. 
Simulation results and conclusions are presented in the end of this paper. 
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OBJECT BASED MODEL 

In this paper, we addresses the problem of detecting high velocity structures in low velocity 
background with sparse data. The relation between the travel time and velocity model is 
illustrated in Fig. [TJ where the forward step predict the travel time for a given velocity 
model. In the inverse step, a velocity model is estimated based on the measured travel time 
data. In this section, we will introduce the object based approach to represent the high 
contrast velocity model. 



Forward step 




Inverse step 

.png 



Figure 1: The forward-inverse step of the velocity model and travel time 

To represent the high velocity structure with the object based model, we need to define 
the type of object that is appropriate for our application. In rendering (computer graphics), 
it has been well studied that we can approximate an arbitrary structure in any detail with 



few fundamental objects (Agarwal and Suri, 1998), for example, triangle or ellipse. In this 



chapter, we use convex objects as our fundamental objects, which can be selected based on 
user's choice. For example, if we only use one kind of convex object (ellipse) as the fun- 
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damental object, the zth object wih be parameterized by the scale (size of object), center 
location and the orientation of the corresponding ellipse, leading to a vector of parameters 
— {<5i, c^, where Si^Ci^il^i denote the scale, center location and the orientation re- 
spectively (see Fig. [2]^b)). We denote \0i\ the area inside the zth object. This formulation 
can be easily extended to multiple convex shapes, for example, the object can be either a 
triangle or an ellipse with one more parameter to identify the type of shape. 

5 - ; 
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15 - ; ; ; ;■■ 

20 - ; ; 

25 - ; ; ; 

30 - ; i ;■■ 

5 10 15 20 25 30 35 40 

(a) (b) 




Figure 2: (a) Grid based model: the HVS structure is approximated with cells and assigns 
the same velocity inside each cell, (b) Object based model with ellipse as the fundamental 
object: the structure is approximated with objects and assigns the same velocity inside each 
ellipse. 

Next we define the velocity model in terms of the objects and the background. Let 
{01,02, • • • ^^n} be the set of N high velocity objects, we denote Vh the velocity in the 
homogeneous background, and Vi the velocity inside object Oi. Therefore, the velocity 
distribution can be represented by the objects as: 

{Vi, if (x,y) e\Oi\, 
(1) 
Vh, otherwise. 
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Obviously the spatial velocity distribution is determined by the objects. We use the notation 
V{Oi^ . . . ^On) to indicate the velocity distribution when we have N objects with parameters 
01, . . . , 07V in the model. 

We use d{Oi, Oj) to represent the distance between two objects, which is defined as: 

(i(0^, 0o) = min ll/x — i/IU, /x E |0J, ly^lOjl (2) 

and the corresponding path is denoted by P(0i,0j). The same notation can be used for 
the distance between point and object, or point to point, i.e., (i(a,/3) and P(a,/3) are the 
distance and path between points a and /3, etc. 



FORWARD STEP 

In this section, we introduce the forward step: Given the velocity model as input, the 
forward model calculates the corresponding travel time between arbitrary transmitters and 



receivers. We use the mathematical formulation proposed in (Berryman, 1991a) to calculate 
the travel time, where the travel path is defined as the direction of wavefront propagation. 
Based on Fermat's principle, the actual travel path is the one with minimum time cost from 
all possible travel paths connecting the transmitter and the receiver. Let V{x,y) represent 
the velocity at position The time cost of an arbitrary path P connecting two points 



{a,/3} based on velocity model V is defined by the path integral: 



p 



r 



(y, a,^) = / dl^, with Pstart OL, Pend f^- 

Jp V{x,y) 



(3) 



We define the travel path, P*, as the path with minimum time cost r*. Therefore, we can 
define the travel time r* between two points {a,/3} as: 

r*(y,a,/3) = minr^(y,a,/3) = min / Ty^^df, (4) 
P P Jp y{.x,y) 
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the path P* will be: 



P*{V,oc,/3) 



argmin r^(V,(x,/3). 

p 



(5) 



Finding the analytical solution for the travel time r* and travel path P* is a classical 



problem in calculus of variations (Courant and Hilbert, 1962). Conventional methods, in- 



cluding the shotgun ray-tracing method (Berryman, 1991a) or the level set method (Sethian 



1996), tend to be all very computationally expensive. In what follows, we will show that 



simple solutions exist if we consider the case of a homogeneous velocity background media 
containing high velocity convex objects. 



Fast travel time/path finding 

For the object based velocity model proposed in the previous discussion, we now prove 
that a travel path finding algorithm can be built using induction on the number of objects 
included in the model. If there are N objects in the model and their velocities are different 
from each other, the computational complexity will be 0(N\). For now, in order to reduce 
the complexity, we assume that all objects have the same velocity and consider the high 
contrast velocity case, i.e., v = vi = • • • = vn and velocity ratio v/v^ oc. With the 
high velocity contrast assumption, we can ignore the time spent passing through an object, 
which provides a way to find the fastest travel path by considering the path between objects 
recursively, thus greatly reducing the computational complexity. 

We build our path-finding algorithm by induction. We start by assuming there is only 
one object in the velocity model, V{Oi)^ then the travel path between a transmitter a and 
a receiver /3 will be either the direct path connecting transmitter-receiver points a and /3 
or the one passing through the object, whichever is faster (see Fig. [3]^ a)). Thus, the travel 
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time function for any two points cx and /3 with a single object 6i in the model is the fastest 
of two possible travel paths: 

l/vh-d{a,P) 

T* (V (Gi), a, 13) = mini (6) 

[ 1/vh ■ ( d{a, 9i) + d{ei,P) ) + 1/vi ■ d(7, C) 

where 7, ^ are the closest points to a, f3 in \6i\, which are defined as: 

7 = argmin ||q; — 7II, ^ = argmin ||/3 — 7,Cg|^i|- (7) 

Because of the homogeneous background and convex object assumption, the travel path 
passing through the object is the combination of the shortest distance from a to |0i|, a 
path inside \6i\ and the shortest distance from \6i\ to /3, which are all straight lines and 
denoted P(a,7), P(7, ^) and P('y,/3) respectively. This travel path can be viewed as 
OL — Oi — j3 which describes the order of objects passing through. From the high velocity 
object assumption, the travel time inside \0i \ is negligible, thus, the travel time function in 
the one-object case becomes: 

l/vh ■ d{a,(3) 



T* {V{ei), a, f3) ^ mm { 



(8) 

l/vh-{d{a,ei) + d{ei,l3) ) 



Next, we consider the two-object case, N — 2. Now the fastest travel path will be either 
the path using less than two objects, or the path passing through both objects. From ([s]) we 
know how to compute the travel path/time for the one-object cases, which are ol — Oi — (3 
and a. — O2 — /3. Thus, all we need to do now is to compute two new paths which pass 
through both \0i\ and \02\^ and compare them to the previous results. 

We notice that the path oc — Oi — O2 — l3 includes the shortest paths from a to object 
01, from object Oi to object 02, and from object O2 to /3. Similar to ([s]), the corresponding 
travel time will be l/v^ • ( d{oc, ^i) + d(0i, ^2) + ^(^2, /3) ). But Oi) has been calculated 
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Transmitter 
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Figure 3: The travel path will be the faster one between the direct path or the one through 
high velocity object 

in the previous step when finding the path a — 0i — /3, and likewise d{02,l3) is also available. 
Thus d{0i,02) is the only new quantity to be calculated. Now it is clear that the travel 
path a. — Oi — O2 — f3 can be found efficiently based on previous results. Similarly for 
cx — O2 — Oi — f3, only (i(02, ^1) would be needed, but (i(02, ^1) = ^(^i, ^2) and thus no new 
distance needs to be computed. 

Now the forward model for calculating the travel time given two objects in velocity 
model T*(y(0i, ^2), yS) can be simplified as 



d(a,0i) + d(0i,^) 



T*(y(0i,02),a,^) = lM-min<^ 



d(a,6>2) + d(6>2,/3) 



(9) 



rf(a,0i) + rf(0i,6>2) + d(6>2,^) 



d{cx,e2) + d{ei,e2) + d{ei,f3) 
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Note that compared to the single object case in ([s]), in ([9]) we only need to compute 
three new terms (i(a,02), ^(^1,^2) and d{02,f3). By induction, it follows that when 
we add the A^th object On the new terms that need to be computed will be d^oc^Ojsi)^ 
d{6i^ Ojsf)^ . . . ^ d{6]s[-i^ On) and d(0N^(3). That is, the number of new distances to be com- 
puted increases linearly with the number of objects. Therefore, the overall computational 
complexity of path tracking becomes 0{N'^). 



Dijkstra path finding 

Note that in Q, as well as in successive cases with additional objects, N > 2, the optimal 
result is obtained by comparing different combinations of pairwise distance (between points 
and/or objects). Based on this fact, we can convert the path tracking problem into a 
shortest distance problem on a graph. 

To do so, we can define a graph where the transmitter and receiver act as source and 
destination vertices and the objects are intermediate vertices. The weight of edges between 
the vertices is the distance between the corresponding objects defined in ^ (sources, des- 
tinations or objects, see Figure [4]). Thus, this undirected graph G = {v^e) will be fully 
connected with nonnegative weights and the shortest path from source to destination can 
be found by running the Dijkstra algorithm (Dijkstra, 1959| ) (see Algorithm [T]) . 



For example, in Fig. |4(b)| we can construct the corresponding fully connected graph for 
velocity model V{0i^02) with 4 nodes as shown in Fig. ^h). We show the update of the 
path tracking algorithm from the source a to each node for a numerical example with the 
same topology as the example of Fig. [ij^b) in Fig. [sj 
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(a) (b) 

Figure 4: Graph representation of path tracking, (a) The distance between objects (b) 
Graph and distance metric 

Relationship between object size and travel time 

For high velocity convex objects, we note that there is a monotonicity property between 
the travel time function and the size of the objects. This will play an important role in the 
inverse step to be presented later. To see this, we assume two high velocity convex objects, 
{^1,^1}, where \0i\ is a dilation of \0i\. In other words, Oi is a convex object with larger 
size but with the same overall shape as 61, thus, \0i\ is a subset of \0\\, \6i\ C \0[\ (as 
shown in Fig. [g]). 

From previous discussion, the travel path depends on the distance between the objects. 
If one object expands, the distance between this object and others must be shorter. Thus, 
all the edge weights (distance metric) in the graph will be smaller than or equal to the 
weights before expansion, thus, the travel time must be faster. This leads to the following 
lemma: 
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(a) 



(b) 




(3)..L.<4) 




(d) 



Figure 5: The example for Dijkstra algorithm. Note in (b) the dist[02] is dist[a]+e(a, O2) 
5, and in (c) after we add Oi the dist[02] is updated to dist[0i] + e{Oi, O2) = 4 
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Figure 6: Change of the travel path with respect to the object size 

Lemma 1. Consider two high velocity objects, {Oi^Oi}, where \6i\ is a dilation of \6i\, 
\0i\ C \6i\. With the other N — 1 objects fixed, consider two alternative velocity mod- 
els V{Oi, . . ^On) and V{6i, 62, • . ,0n) using 61 and 61 respectively. For the travel 
time function, we have T*(y(0i, . . . , On), a, /3) > T*(y(0]^, . . . , On), a, /3) Va, /3. Thus, the 
travel time is monotonically non-increasing with respect to the size of high velocity object, j^] 

Proof. The proof is straight forward, given that weights in the graph are smaller or equal 
with the same topology, the shortest path will be shorter. □ 

INVERSE STEP 

In the forward step we just presented, we can predict the travel time when the velocity 

model is given. Then, in the inverse step the goal is to estimate the velocity model when the 

travel time data is observed. With limited travel time measured data, this inverse problem 
*This property still holds for any even \6i \ G \0i\, even \6i\ and \6i \ have different shapes. 
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becomes ill-posed and the solution may not be unique. In this case, searching one single 
most possible model provides little information for the solutions of this inverse problem. 
Instead, we formulate it as a statistical inference problem and estimate the probability 



distribution of the velocity model in the parameter space (Tarantola, 1987). We start this 
section by introducing some notations. 



Notations 

The input data is obtained by measuring the travel time between the set of transmitters 
A = {ai, . . . , cxtx} and receivers B = {/3i^ . . . ^/S^x}- We denote the measured travel time 
for all transmitter-receiver pairs {cxi^/3j) as a vector t = {ti, . . . where n — tx - rx. 
Assuming that there are at most N objects in the velocity model, we can cascade all object 
parameters and define a vector of model parameters, = . . . , 9n}^ thus the velocity 
model V{Oi, . . . ,0n) can be represented simply by 1/(0). Then we define the travel time 
function from the forward model, T(@,A,B), as a vector function representing the travel 
time between each pair of transmitters and receivers based on the velocity model with 
parameter 0: 

T(0, A, B) = (ri(0, yl, S), . . . , Tn{@, A, B)), (10) 

where 

Ti{@,A,B) = t*{V{9i,...,9n),cxuI3^) 

\ (11) 

Tn{@,AB) = T*{V{ei,...,eN),CXtx,l3r.)- 

We then define an error function as a quadratic data-fitting error between the travel time 
predicted from the forward model and the measured travel time: 

E{@)^\\t-T{@,A,B)f. (12) 
19 



We use the Bayesian approach (iTarantolal 119871), which leads us to update the behef for 



different models after accounting for the observations. The likelihood function L{@), which 
is modeled as a Gaussian uncertainty of the error, measures the confidence on different 
models: 

L(0) = fc.e-^(®), (13) 

where A: is a normalization constant. From Bayes' rule, the posterior probability density 
function (PDF) cr(0) is proportional to the prior probability distribution p{@) multiplied 
by the likelihood function L(@). The prior probability p(0), which may come from pre- 
vious experience, can provide useful information to select possible models after the data is 
observed: 

a{@)^k-p{@)L{@), (14) 

where k is again a normalization constant. For the rest of this paper, we assume a uni- 
form distribution for the prior probability p(0), so that the posterior PDF is equal to the 
likelihood function. Estimate the posterior PDF is equivalent to estimating the error func- 
tion, where the most possible models correspond to the global minima in the error function 
E{@). In the next section, we will introduce our algorithm to estimate the error function. 



Proposed algorithm 

Because the travel time function is nonlinear, the error function £"(0) will have a complex 
and multi-modal shape. We can "sample" the error function at any point by calculating 
the value of travel time function from forward model for the corresponding parameters. 
However, even with the fast path tracking approach of Algorithm [T} calculating the travel 
time is still very computationally expensive. 
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A naive approach to estimate the error function would be running a brute-force uniform 
samphng in the whole parameter space. For example, consider the 2D case where each 
object needs one parameter for the size, two parameters for the center location and one for 
the rotation angle. If we choose 3 objects in the velocity model, there will be 12 parameters, 
so that if we uniformly sample 10 possible values for each parameter, this will require a total 
of 10^^ samples, which clearly makes this uniform sampling approach impractical. 



We note that in (13) the high probability models correspond to the regions with low 
error in the parameter space. Thus, when sampling the error function, we would like to 
have more sample points in these low error regions. However, the error function is multi- 
modal and the gradient based method (steepest descend) can only search and sample near 
the closest local minima. Thus, in order to search and sample in the whole parameter 
space, we choose a random walk sampling scheme. While this is a popular approach to find 
multiple minima, the main drawback is that its computational time could be very high. 



especially when the dimension of the parameter space is high (MacKay, 2003). Thus, when 
we consider a velocity model with many objects, the dimension of the model parameter 
space grows with the number of objects which implies an exponential growth in the number 
of samples. To overcome this problem, we propose an accelerated sampling algorithm to 
speed up the random walk sampling and achieve denser sampling in the low error regions. 

To further accelerate the sampling, we make use of known properties of the error func- 
tion. Specifically, in the second part of our algorithm we use the monotonicity property 
in Lemma [l] which implies that each error function (corresponding to one measured travel 
time) is a unimodal function with respect to changes in the size of one object. This allows 
us to re-sample along the dimension corresponding to the object size (with all other pa- 
rameters fixed). These re-sampling location are chosen by running golden section search on 
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each error function (which is unimodal). Compared to other samphng methods, our method 
significantly reduces the computational time and provides a sufficiently good approximation 
of the error function. The two parts of our algorithm are described in the next section. 



Accelerated random walk sampling 

In the first part of our algorithm, we want to sample the error function and emphasize 
the sampling in the low error regions. We modify the "Hamiltonian Monte Carlo" (HMC) 



method (Duane et al., 1987), which is a Metropolis method but includes the gradient infor- 
mation to reduce the random walk behavior. HMC uses the dynamical system concept to 
draw samples by simulating a particle movement on the surface of the error function. We 
introduce a new "friction" term, which is closely related to the cooling schedule in simulated 



annealing (Kirkpatrick, 1984), to draw more samples near the local minima. 



In HMC, we define a dynamical system where the model parameter is augmented 
by a momentum variable p, where p, have the same size. The total energy i7(0,p) of 
the dynamical system is defined as the sum of "kinetic energy" and the "potential energy" , 
where the "potential energy" is equal to the error function £"(0) and the "kinetic energy" 
is given by K{p) = ||p|p/2, i.e, i?(0, p) = £"(0) + K{p). The changes in and p will be 
determined by the following equations: 

e = p, (15) 

P = -^-P- (16) 
To sample the error function, we simulate and record the states of particle movement. 



With a randomized momentum p, we solve the Hamiltonian dynamics in (15) during a 
simulated time of duration t. In the dynamical movement, the momentum variable p 
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determines where the parameter goes, and the gradient of potential function determines 
the change of momentum p. The friction term ep decides the loss of total system energy. 

Starting with the initial value ©o, we define the simulation time t and steps At, then 



use (15) to identify the successive steps in the walk through parameter space. During the 
simulated time t, we can record the state variables [©(At), p(At)], [©(2At), p(2At)], . . . , 
[©(t),p(t)] which describes the variable movement on the error function. Then we take all 
the ©(At), . . . , ©(t) as new sampling points for the error function. 

In Fig[7|we show an example of the particle movement. The initial momentum drives the 
particle to the high potential region, then it falls back because the momentum is changed by 
the gradient of the potential function. And the total energy decreases during the simulation, 
which causes the particle to settle down in one local minima. 

Then we decide whether to accept the last sampling point ©(t) as a new starting point 



for the next round of simulation by the Metropolis rule (Metropolis et al., 1953). This 
acceptance rule is based on the change of the error function, A£'(©) = E(@(t)) — £"(©(0)), 
where ©(0) is the current starting point. If A£'(©) < 0, the new sample reaches a lower 
error state and we always accept it as a new starting point. Otherwise, we draw a random 
number r between [0, 1] and accept it if exp(— A£^(©)) > r. We iteratively simulate this 
dynamic system L times, and the detail algorithm is shown in Algorithm |2j 

This sampling scheme provides some useful properties for exploring the parameter space 
in our inverse problem: 

• The random momentum provides the ability to jump out of current local minima, and 
makes it possible to travel through all parameter space and sample in multiple low 
error regions. 
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Figure 7: The simulated state movement. Note the initial momentum is toward left, driving 
to the high error region. Then it falls into low error regions, and we can see the total energy 
is decreasing through the simulation. Most of samples are near local minimum. 
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• The additional friction term decreases the total system energy at each simulation 
proceeds, which "cools down" the system and allows the state to stay near a local 
minimum. This property leads to more samples in the low error region. 

• The exploration speed is linearly related to the number of iterations, instead of being 
related to its square root, as in the typical random walk. This makes our sampling 
scheme much more efficient in a high- dimensional parameter space. 

Re-sampling by the monotonicity property 

In the first part of our algorithm, an accelerated random walk leads to sample points 
concentrated in the low error regions. However, if we want to understand the structure of 
the error function in the whole parameter space, it will not be sufficient if we only sample 
in the low error regions. In the second part our algorithm, we want to re-sample the error 
function and build a linear approximation for it. 

Because the shape of error function is complicated and multi-modal, without any prior 
information there is no easy way to efficiently choose "good" sampling locations to approx- 
imate it. One possible approach would be using the samples from HMC step as starting 
points to run a pure random walk to explore and re-sample the error function. However, 
the exploration speed in a pure random walk would be roughly equal to the Nth root of the 
number of samples, where N is the number of parameters. This would be too slow and we 
would require exponential samples to cover the whole parameter space. We now show how 
to use some properties of the error function in order to choose the re-sampling locations 
efficiently. 

We note that the error function is defined by the sum of mismatches corresponding to 
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each measured travel time: 

E{@) = ||t-T(0,^,e)f (17) 
= pi - ri(0, A, B)f + \\t2 - r2(0, A,B)f + ... (18) 
= Ei{@) + E2{&) + ..., (19) 

where each separate error function has a quadratic form Ei{Q) = \\ti — T^(0, ^, /3)|p. 
Consider the change in the travel time function Ti{@^A^B) with respect to the size sj 
of object j, while all other parameters are left unchanged. By Lemma [l| the travel time 
between two arbitrary points will always be non-increasing as the size of a high velocity 
object increases. Thus, given its quadratic form the change of each individual error function, 
£■^(0), with respect to a single object size will be a weakly unimodal function (see Figjs]). 




Size of object Size of object 

(a) (b) 



Figure 8: Change of travel time and error function with respect to the size of one high 
velocity object, (a) The travel time function Ti{Q,A^B) (b) The error function Ei{Q = 
\\ti — Ti{@,A,B)\\'^. The travel time is monotonically non-increasing, therefore the error 
function is weakly unimodal. 
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We make use of this property to choose the re-samphng locations in the parameter space. 
The basic idea is to re-sample only along the dimensions corresponding to size parameters, 
si, 52, . . . . We choose one size parameter each time and use the above property. Because 
Ei{@) is unimodal with respect to the change of 5j, it can be well represented by linear 
interpolation if we sample more frequently in large curvature regions. Thus, we select one 
axis Sj at a time corresponding to the size of object j and use the golden section search 



(Kiefer, 1953) to perform selection on re-sampling locations along that axis, while leaving 
parameters along the other dimensions unchanged. We use these probing points as re- 
sampling locations, with more locations chosen in the large curvature regions (minima). 
Thus, compared to the random or uniform sampling, the golden section search sampling is 
more efficient because it puts few samples in the almost constant regions and focuses on 
the large curvature regions in the parameter space. 

To illustrate our approach, consider an example where we have two objects in the velocity 
model = {01, 02 } and two measured data points t = {^1,^2}, the error function is defined 
as £"(0) = £"1(0) + £2(0)- Given an initial sample 0(0), the parameters can be re- 
ordered as 0(0) = {51(0), 52(0), 0(0)}. To find the re-sampling locations on the si axis, 
we run the golden section search k times along the si dimension for each of Ei(@) and 
£2(0). The re-sampling locations wih be {5i(l), 52(0), 0(0)}, . . . , {5i(A:), 52(0), 0(0)} and 
{si{k + 1), 52(0), 0(0)}, {5i(2A:), 52(0), 0(0)} where the first k points are chosen by 
running the golden section search on Ei and the next k points are on E2. Likewise, we run 
the same procedure to choose the re-sampling locations on the 52 axis. 

Note that we choose 2k re-sampling locations for the error function on the axis si based 
on El and E2 separately. The reason is that the error function is the sum of each separate 
error function £"(0) = £i(0) + £2(0)5 which may not be unimodal (sum of unimodal 
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functions is not necessary unimodal, see Fig. |9]). Thus, running the golden section search 
directly on E(@) may not give us good sampling locations. Since each separate function is 
unimodal and we know how to efficiently sample it, we divide the sampling "budget" and 
choose the sampling locations based on each separate function. Because the error function 
is the sum of each separate error function, the large curvature regions for the error function 
should belong to the union of large curvature regions of each separate error function and 
our approach selects more samples in the regions where one of the two error functions has 



large curvature. In Fig. [TO] we show an example of how this approach works better than 
choosing sampling locations directly on the total function. 





(a) 



(b) 



Figure 9: Sum of two unimodal functions. Note the sum of two unimodal functions is not 
necessary a unimodal function. 

SIMULATION RESULTS 

Following our previous assumptions, in our simulations we use a velocity ratio v/v^ ^ oo 
and choose "line" as the pre-defined convex geometrical object to model the HVS. The model 
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(a) (b) 

Figure 10: Re-sampling comparison, (a) Re-sampling locations chosen by running the 
golden section search on f(x) (b) Locations chosen separately on and f2{x). The red 

"squares" are sampling locations chosen from f2{x), and black "dots" are from fi{x). Note 
that choosing from separate functions gives much better approximation because each one 
is unimodal. 
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parameters Oj = {5j,Cj,'0j} will be the length, line center and angle for the line object. 
The scenario we choose for our experiments is motivated by the problem of modeling a 
fractured reservoir, where the fractures are usually represented by lines in 2D or planes in 
3D. 

To visualize the object based model in spatial domain, after we sample {©i, . . . , Qn} 

N points in the parameter space, we define a mapping function /(©) which maps the object 

parameters into object shapes in spatial domain. Then we calculate the average: 

N 

Mf= a{&)f{&)d& ^ 5^a(0,)/(0i). (20) 

Because /(©) represents the high velocity structure in spatial domain, A4f can be viewed 
as the "appearance probability map" of the high velocity in different region. 

In simulations, we need to choose the random momentum p and friction e • p. We 
draw the momentum from a normal distribution A^(0,p), where p is equal to the 10% of 
the maximum possible value of 0. The friction coefficient e is chosen to be 1/t in our 
experiments . 

In Experiment 1, we illustrate the sampling process. Assuming there is only one 
measured data point (one travel time between transmitter-receiver pair), and the center of 
"line" object is fixed. Thus, in this case we only have two parameters {s^tp} which represent 
the length and angle of the line object. We show the geometry of the line object, sensor 
locations and the ground truth PDF (which is multi-modal because only one measured 



data point) in Fig. 11, Now we apply our algorithm to sample and estimate the PDF: 
The samples from the HMC step are listed as blue circles, which are concentrated in high 
probability regions. The results of randomized re-sampling and golden section search are 



shown in Fig. [T2| where the randomized re-sampling only explores a very narrow region 
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near the starting points and golden section search re-samples along the whole s axis. 



The estimated PDF from the two different approaches are shown in Fig. [T3| which 
shows that our proposed method has much better estimated PDF than the result from 
randomized re-sampling. Because the randomized re-sampling only explores a small region 
near the initial samples, most of the structure of the PDF remains unknown. Our approach 
choose the re-sampling locations which lead to better approximation of PDF structure. 
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(c) 



Figure 11: The (a) geometry of sensor settings and (b) (c) the PDF for experiment 1. 

In Experiment 2, we increase the number of measured data (2 transmitters and 2 
receivers, total 4 measured travel time between transmitter-receiver pairs). The PDF and 
our estimation are shown in Fig. [T4j Comparing to the Experiment 1, the change of PDF 
is much sharper which implies the model uncertainty is less if we have more data. For 
example, if we define the acceptable model as £"(0) < 5, then the corresponding region in 
the parameter space is much narrow for the large data case. 



We show the ground truth model and the estimated probability map in Fig. 15 The 
result shows that we have high probability areas near the ground truth. Note the center has 
highest probability - the reason is that in this simple example we fix the center for all HVS 
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(b) 



Figure 12: The HMC sampling and randomized re-sampling in experiment 1. (a) The blue 
"dots" are samples from HMC. (b) The red "x" are re-samples from random walk. Note it 
can only explores a small region in parameter space, (c) The red "x" are re-samples from 
golden section search. The re-sampling locations have the same il; but cover the whole s 
axis. 




10 15 20 25 30 35 40 45 50 




(a) 



(b) 



Figure 13: The estimated PDF by (a) Randomized re-sampling (b) Golden section search 
re-sampling. Note the randomized re-sampling only explores a small region of parameter 
space and most of PDF is unknown. 
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models. Thus, when we calculate the appearance probability, all models have a common 
center point and we will have the highest probability near the center area. 




Figure 14: The PDF and our estimation in experiment 2. Note comparing to experiment 
1, the PDF has sharper changes which implies higher model resolution. 

In Experiment 3, the transmitters and receivers are placed on the boundary and the 
HVS is near the center lower area. In the inverse step, we use one line object in the HVS 
model and all object parameters (center position, length, angle) can be changed freely, so 
that the parameter space has 4 dimensions. 

To quantify the reconstructed model error, we define the error metric as the difference 
between each sampled HVS model and the ground truth, weighted by the probability of 
each estimated HVS model. We use the Hausdorff distance to characterize the difference of 
two HVS models. The Hausdorff distance for two objects X, Y is defined by: 

(iif (X, Y) = max < sup inf (i(a,/3), sup inf a) > , (21) 

where (i(a,/3) is the distance function of Q. 
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Figure 15: (a) The ground truth and (b) the appearance probabihty map. It shows several 
different models closed to ground truth all have high probability. 

In our result, we map the high velocity model into 2D spatial domain, then discretize the 
result into a set of points. The Hausdorff distance can be viewed as the maximum distance 



among all points in a set to the nearest point in the other set (Rote, 1991) and is widely 



used in computer vision to measure the difference between 3D curves or binary images (Hut 



tenlocher et al., 1993). It also has an interesting property in that when (i/f(X, Y) = then 
X and Y have the same closure. In our case, because we use convex object models, if the 
Hausdorff distance is zero it implies that the two objects are equal X = Y. From previous 
discussion, for a given model parameter 0, /(©) is the function that represents the shape 
of the corresponding high velocity structure in spatial domain. Then, the error metric with 
respect to the true HVS ©truth can be written as 



<f(0) = J dH{f{@)J{@truth))-<j{@)d@ 

N 

~ J]d/f(/(©i),/(etruth))-<7(0i), 



(22) 

(23) 



1=1 
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where 0^ are sampled models. 



We list the error metric for 5*5, 10*10, 25*25, 50*50 transmitter-receiver pairs (see Fig. 
16). Our results show that the error metric decreases as the number of sensors increases. 
We also list the result with grid-based model and linearized reconstruction algorithm for 
comparison. Due to the travel paths only covering very few cells, most of the cell's velocity 
remain unknown and the linearized reconstruction algorithm will assign the high velocity 
to cells which give the most significant changes for the error function, which usually are the 



cells closed to sensors, (see Fig. 17) 
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Figure 16: (a) The ground truth and (b) the error metric for experiment 3. Note that the 
error metric decreases when the number of sensor increases. 

In Experiment 4, we use the same sensor constellation as in Experiment 3 with a more 
complicated HVS. For the inversion, we use 3 line objects and the result shows some diffi- 



culties to resolve the vertical structure of HVS. In Fig. [191 and [20} shows that increasing 
the number of measured data points does not increase the vertical resolution. This is an 
inherent limitation of travel time tomography, which comes from the relative location of the 
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Figure 17: (a) The grid-based result with 25*25 transmitters/receivers and (b) the appear- 
ance probabihty map for Experiment 3. Note the grid-based model fails to recover the HVS 
location. 

sensors and the HVS. If the "line" HVS is orthogonal to the travel path, it will not affect 
the travel time at all. In this case most of the transmitter-receiver pairs are in the hori- 
zontal direction, so that the vertical resolution is very limited. This is the reason we have 



"phantoms" in the vertical direction and the error metric (see Fig. 18) does not decrease 
when we increase the number of measured data points. 

From the above experiments, we show that our algorithm can estimate the PDF and 
recover possible models for high contrast travel time tomography with sparse data. Our 
algorithm is robust to noise because we estimate the PDF from all measured data. For 
example, if one measured travel time is affected by shot noise, the corresponding error 
function Ei will be significantly changed. But the error function is defined as the sum of all 



individual error functions (see equation ( 12 )) , thus, the effect of shot noise will be "averaged" 
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Figure 18: (a) Ground truth and (b) error metric for experiment 4. Note the error metric 
does not decrease monotonicahy due to the vertical phantom. 
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Figure 19: (a) Grid-based result and (b) appearance probability map for 10 * 10 sensors in 
experiment 4 
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(b) 



Figure 20: (a) Grid-based result (b) Appearance probability map for 25 * 25 sensors in 
experiment 4 

out by other good measure data points which provides the robustness for estimation. 

CONCLUSION 

The main purpose of this paper is to propose a new approach for the reconstruction of high 
contrast discrete velocity models in travel time tomography. To image the high contrast 
media, we take advantage of the high velocity structures and prove that the travel trajectory 
is piecewise linear when we consider convex object model. We show that our travel path 
finding method is able to compute the corresponding travel time much faster than the 
conventional ray-tracing method because we consider the number of objects as "nodes" 
instead of number of cells in grid model. Thus, the complexity scales with the number of 
objects, instead of the (much larger) number of cells in a grid. 

A model based approach for high velocity structures (HVS) is studied and the error 
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function is defined by the misfit between the predicted travel time based on current model 
and the measurements. We develop a reconstruction algorithm to efficiently sample the error 
function in the model parameter space. After we map the corresponding model parameters 
back to the object shape in spatial domain, we can obtain the HVS appearance probability 
in different areas. In our simulations, we show how our algorithm samples and approximates 
the error function, finding set of possible HVS models. The results also show our algorithm 
has better reconstructions compared to the typical grid-based model. 

Future work will be to explore how to refine the idea of object based model to achieve 
an efficient representation of the HVS structure. For example, we are looking at how to 
define or adaptively change the shape of objects to have a sparse representation. And the 
optimum number of objects for the model is another interesting question. Increasing the 
number of objects will provide a better detail representation of structure, but also increase 
the dimension of parameter space. We plan to explore the trade-off between the number of 
objects and computational complexity of our randomized sampling algorithm. 
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Algorithm 1 Dijkstra algorithm for path tracking 



for V e G do 

dist[7;] = oo 
previous ['L'] = 
dist [source] = 
Q := G G 
end for 



> Initiahzation 
t> Unknown distance from source to v 
> Previous node in optimal path from source 
> Distance from source to source 
> Put all nodes in Q to be scanned 



> The main update loop 

> Start node in first case 

> remove u from Q 



while Q / do 

u := V e Q with minimum distf^;] 
Q = Q\u 
if dist[t^] = oo then 

break ; > all remaining vertexes are inaccessible from source 

else 



for V neighbor v of u do 
alt = dist[t^] + e{u^ v) 
if alt < dist['L'] then 
dist[7;] = alt 
previous ['L'] = u 
end if 



> where v is still in Q 



> update the distance for v 



end for 
end if 
end while 
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Algorithm 2 Modified Hamiltonian Monte Carlo method 



©(0) = @^n^t 

for I = 1 : L do 

9 = v£;(0(o)) 

p ^ randn( size(0(O)) ) 

®new — ®(0)5 Snew — S 

for tsim = : At : t do 

p = p- S ' gnew/"^ 
®new — ®new + ' P 

Sample list ^ ®new 

p = p- 5 ' - e • p 

end for 

^new — E{®new) 
AE = Enew — E 

if AE <0 then 



> Initialize 
> loop L times 
> set gradient using initial parameter 
> set error function value 
> initialize momentum 

> Use "leapfrog" steps to simulate the dynamics 



t> Record the state samples 
> Update the gradient and momentum 



t> Find the final value of error function 



> Use Metropolis rule to decide the new starting point 



0(0) = ©(*) 



> Accept the new starting point 



else 



if exp{—AE) > rand() then 

0(0) = @{t) t> Accept the new starting point if exp(— A£^) < a random 



number 

end if 
end if 
end for 
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